function [A, b] = process_boundary_2d(A, b, mesh, pde)
ref_gauss_points_1d = [-0.9681602395076260898355762, -0.8360311073266357942994298, -0.6133714327005903973087020, -0.3242534234038089290385380, 0, 0.3242534234038089290385380, 0.6133714327005903973087020, 0.8360311073266357942994298, 0.9681602395076260898355762];
ref_gauss_weights_1d = [0.0812743883615744119718922, 0.1806481606948574040584720, 0.2606106964029354623187429, 0.3123470770400028400686304, 0.3302393550012597631645251, 0.3123470770400028400686304, 0.2606106964029354623187429, 0.1806481606948574040584720, 0.0812743883615744119718922];
psi = reference_basis_function_2d(mesh.basis_type);
%% Robin
matrix_A = sparse(mesh.Nb,mesh.Nb);
vector_b = zeros(mesh.Nb,1);
for k = 1:size(pde.bdary,2)
    if pde.bdary(1,k) == 2
        bde = boundaryedges(mesh, k);
        switch k
            case 1
                cq = @(x,y) pde.cor(x,y).*(-pde.dudy(x,y) + pde.r(x,y).*pde.u(x,y));
                cr = @(x,y) pde.cor(x,y).*pde.r(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),1);
                    p2 = mesh.node(bde(3,n),1);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cq(gauss_points_1d(j),mesh.bottom)*psi{i}(x_loc(gauss_points_1d(j),mesh.bottom),y_loc(gauss_points_1d(j),mesh.bottom));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                        for j = 1:mesh.Nlb
                            r = 0;
                            for beta = 1:9
                                r = r + gauss_weights_1d(beta)*cr(gauss_points_1d(beta),mesh.bottom)*psi{i}(x_loc(gauss_points_1d(beta),mesh.bottom),y_loc(gauss_points_1d(beta),mesh.bottom))*psi{j}(x_loc(gauss_points_1d(beta),mesh.bottom),y_loc(gauss_points_1d(beta),mesh.bottom));
                            end
                            matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) = matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) + r;
                        end
                    end
                end
            case 2
                cq = @(x,y) pde.cor(x,y).*(pde.dudx(x,y) + pde.r(x,y).*pde.u(x,y));
                cr = @(x,y) pde.cor(x,y).*pde.r(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),2);
                    p2 = mesh.node(bde(3,n),2);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cq(mesh.right,gauss_points_1d(j))*psi{i}(x_loc(mesh.right,gauss_points_1d(j)),y_loc(mesh.right,gauss_points_1d(j)));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                        for j = 1:mesh.Nlb
                            r = 0;
                            for beta = 1:9
                                r = r + gauss_weights_1d(beta)*cr(mesh.right,gauss_points_1d(beta))*psi{i}(x_loc(mesh.right,gauss_points_1d(beta)),y_loc(mesh.right,gauss_points_1d(beta)))*psi{j}(x_loc(mesh.right,gauss_points_1d(beta)),y_loc(mesh.right,gauss_points_1d(beta)));
                            end
                            matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) = matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) + r;
                        end
                    end
                end
            case 3
                cq = @(x,y) pde.cor(x,y).*(pde.dudy(x,y) + pde.r(x,y).*pde.u(x,y));
                cr = @(x,y) pde.cor(x,y).*pde.r(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),1);
                    p2 = mesh.node(bde(3,n),1);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cq(gauss_points_1d(j),mesh.top)*psi{i}(x_loc(gauss_points_1d(j),mesh.top),y_loc(gauss_points_1d(j),mesh.top));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                        for j = 1:mesh.Nlb
                            r = 0;
                            for beta = 1:9
                                r = r + gauss_weights_1d(beta)*cr(gauss_points_1d(beta),mesh.top)*psi{i}(x_loc(gauss_points_1d(beta),mesh.top),y_loc(gauss_points_1d(beta),mesh.top))*psi{j}(x_loc(gauss_points_1d(beta),mesh.top),y_loc(gauss_points_1d(beta),mesh.top));
                            end
                            matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) = matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) + r;
                        end
                    end
                end
            case 4
                cq = @(x,y) pde.cor(x,y).*(-pde.dudx(x,y) + pde.r(x,y).*pde.u(x,y));
                cr = @(x,y) pde.cor(x,y).*pde.r(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),2);
                    p2 = mesh.node(bde(3,n),2);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cq(mesh.left,gauss_points_1d(j))*psi{i}(x_loc(mesh.left,gauss_points_1d(j)),y_loc(mesh.left,gauss_points_1d(j)));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                        for j = 1:mesh.Nlb
                            r = 0;
                            for beta = 1:9
                                r = r + gauss_weights_1d(beta)*cr(mesh.left,gauss_points_1d(beta))*psi{i}(x_loc(mesh.left,gauss_points_1d(beta)),y_loc(mesh.left,gauss_points_1d(beta)))*psi{j}(x_loc(mesh.left,gauss_points_1d(beta)),y_loc(mesh.left,gauss_points_1d(beta)));
                            end
                            matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) = matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) + r;
                        end
                    end
                end
        end
    end
end
A = A + matrix_A;
b = b + vector_b;
%% Neumann
vector_b = zeros(mesh.Nb,1);
for k = 1:size(pde.bdary,2)
    if pde.bdary(1,k) == 1
        bde = boundaryedges(mesh, k);
        switch k
            case 1
                cp = @(x,y) pde.cor(x,y).*(-1).*pde.dudy(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),1);
                    p2 = mesh.node(bde(3,n),1);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cp(gauss_points_1d(j),mesh.bottom)*psi{i}(x_loc(gauss_points_1d(j),mesh.bottom),y_loc(gauss_points_1d(j),mesh.bottom));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                    end
                end
            case 2
                cp = @(x,y) pde.cor(x,y).*pde.dudx(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),2);
                    p2 = mesh.node(bde(3,n),2);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cp(mesh.right,gauss_points_1d(j))*psi{i}(x_loc(mesh.right,gauss_points_1d(j)),y_loc(mesh.right,gauss_points_1d(j)));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                    end
                end
            case 3
                cp = @(x,y) pde.cor(x,y).*pde.dudy(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),1);
                    p2 = mesh.node(bde(3,n),1);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cp(gauss_points_1d(j),mesh.top)*psi{i}(x_loc(gauss_points_1d(j),mesh.top),y_loc(gauss_points_1d(j),mesh.top));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                    end
                end
            case 4
                cp = @(x,y) pde.cor(x,y).*(-1).*pde.dudx(x,y);
                for n = 1:size(bde,2)
                    nk = bde(1,n);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,:);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,n),2);
                    p2 = mesh.node(bde(3,n),2);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for i = 1:mesh.Nlb
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j)*cp(mesh.left,gauss_points_1d(j))*psi{i}(x_loc(mesh.left,gauss_points_1d(j)),y_loc(mesh.left,gauss_points_1d(j)));
                        end
                        vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                    end
                end
        end
    end
end
b = b + vector_b;
%% Dirichlet
for k = 1:size(pde.bdary,2)
    if pde.bdary(1,k) == 0
        bdn = boundarynodes(mesh, k);
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        if pde.err
            b(bdn) = pde.u(mesh.node(bdn,1),mesh.node(bdn,2));
        else
            b(bdn) = pde.D(mesh.node(bdn,1),mesh.node(bdn,2),k);
        end
    end
end
end

%% boundarynodes
function bdn = boundarynodes(mesh, index)
switch index
    case 1
        bdn = 1:mesh.Ny+1:mesh.Np;
        if mesh.basis_type == "P2"
            bdn = [bdn, mesh.Np+mesh.Ny+1:3*mesh.Ny+1:mesh.Nb];
        end
    case 2
        bdn = mesh.Np-mesh.Ny:mesh.Np;
        if mesh.basis_type == "P2"
            bdn = [bdn, mesh.Nb-mesh.Ny+1:mesh.Nb];
        end
    case 3
        bdn = mesh.Np:-(mesh.Ny+1):mesh.Ny+1;
        if mesh.basis_type == "P2"
            bdn = [bdn, mesh.Nb-mesh.Ny:-(3*mesh.Ny+1):mesh.Np+1];
        end
    case 4
        bdn = mesh.Ny+1:-1:1;
        if mesh.basis_type == "P2"
            bdn = [bdn, mesh.Np+mesh.Ny:-1:mesh.Np+1];
        end
    otherwise
        error("Invalid edge number");
end
end

%% boundaryedges
function bde = boundaryedges(mesh, index)
switch index
    case 1
        bde = [1:2*mesh.Ny:mesh.Ne;
            1:mesh.Ny+1:mesh.Np-mesh.Ny-1;
            mesh.Ny+2:mesh.Ny+1:mesh.Np];
    case 2
        bde = [(mesh.Nx-1)*2*mesh.Ny+2:2:mesh.Ne;
            mesh.Np-mesh.Ny:mesh.Np-1;
            mesh.Np-mesh.Ny+1:mesh.Np];
    case 3
        bde = [2*mesh.Ny:2*mesh.Ny:mesh.Ne;
            mesh.Ny+1:mesh.Ny+1:mesh.Np-1;
            2*mesh.Ny+2:mesh.Ny+1:mesh.Np];
    case 4
        bde = [1:2:2*mesh.Ny-1;
            1:mesh.Ny;
            2:mesh.Ny+1];
    otherwise
        error("Invalid edge number");
end
end